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Abstract 

The systems of algebraic equations which arise from spectral 
discretizations of elliptic equations are full and direct solutions of them 
are rarely feasible. Iterative methods are an attractive alternative because 
Fourier transform techniques enable the discrete matrix-vector products to be 
computed with nearly the same efficiency as is possible for corresponding but 
sparse finite difference discretizations. For realistic Dirichlet problems 
preconditioning is essential for acceptable convergence rates. A brief 
description of Chebyshev spectral approximations and spectral multigrid 
methods for elliptic problems is given. A survey of preconditioners for 
Dirichlet problems based on second-order finite difference methods is made. 
New preconditioning techniques based on higher order finite differences and on 
the spectral matrix itself are presented. The preconditioners are analyzed in 
terms of their spectra and numerical examples are presented. 
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1. Introduction 


Spectral methods Involve representing the solution to a problem in terms 
of a truncated series of smooth global functions which are known as trial 
functions. Specifically, these functions are eigenfunctions of a singular 
Sturm-Llouville problem (Gottlieb and Orszag (1977)). This global character 
distinjpiishes spectral methods from finite difference and finite element 
methods. It is also responsible for their superior approximation properties, 
yielding accurate solutions with substantially fewer grid points than are 
required by finite difference methods. 

Test functions are used to minimize the residual that results from the 
substitution of the series expansion of the solution into the differential 
equation. The choice of test functions distinguishes between the spectral 
Galerkin and spectral collocation methods. In the Galerkin approach the test 
functions are the same as the trial functions whereas in the spectral 
collocation method the test functions are shifted Dirac delta functions. For 
many problems, especially nonlinear ones, the spectral collocation method is 
the easiest to implement, and the most efficient as well (Hussalnl, Koprlva, 
Salas and Zang (1983)). The present discussion is confined to the spectral 
collocation method, and all future references to the spectral method will mean 
the spectral collocation method. 

The matrices representing the discrete spectral collocation operator 
corresponding to elliptic problems (Zang, Wong, and Hussaini (1982)) are 
usually full. Even in the constant coefficient linear case, direct inversion 
of these matrices is usually expensive. Iterative schemes are a practical 
necessity. The condition number of the spectral matrices is large, and 
therefore effective preconditioning is necessary. 
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In recent years research in the field of elliptic equations has focused on 
multigrid methods. Basically, the multigrid method is a numerical strategy to 
solve partial differential equations by switching between finer and coarser 
levels of discretization. The characteristic feature of the method is the 
combination of a smoothing step and a coarse grid correction (Brandt (1977), 
Hackbusch (1980)). During the smoothing step the residuals are not 

necessarily decreased but smoothed. In the following correction step the 
discrete solution is improved by means of an auxiliary equation on a coarser 
grid. This results in an iterative method that is usually very fast and 
efficient . 

2. Formulation of the Problem 

We consider the self-adjoint elliptic equation 

^a(x,y,u) + ^b(x,y,u) = f(x,y) (2.1) 

with Dirichlet boundary conditions in the domain [-1,1] x [-1,1]. A proper 
representation of the solution to this Dirichlet problem employs Chebyshev 
polynomials. The details of implementing Chebyshev collocation methods for 
this type of problem have been given by Zang, Wong, and Hussaini (1982) and by 
Hussalni, Salas, and Zang (1983). This discretization of equation (2.1) may 
be written as 

Lgp(V) = F, (2.2) 

where V is 'the vector of unknowns at the collocation points, F is the 

vector of the values of the right-hand side at the collocation points and 
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L is the vector-valued (generally nonlinear) operator representing a 

spectral discretization of the left-hand side. 

For nonlinear problems iterative schemes for (2.2) are a necessity. 
Iterative methods appear to be preferable to direct methods even for constant 
coefficient linear problems since the matrices representing Lgp are full. 
No fast direct methods are known for these systems. Iterative methods are 
appealing because the standard implementation of spectral discretizations 
employs Fast Fourier Transforms which reduces the cost of evaluating the left- 
hand side of (2.2) to 0( N log N) operations, even for nonlinear problems, 
where W is the total number of unknowns. Iterative methods also have a 
clear advantage over direct methods in terms of storage. 

Many iterative schemes can be described within the framework of defect 
correction. Let H be some approximation to the Jacobian J-^ of Lgp(V) and 
let be the latest approximation to V. The simplest iterative scheme is 

the preconditioned Richardson's method 

yn+l ^ ^n _ ^ ^-Ir n^ _ ^2.3) 

where o)^ is a relaxation parameter. 

In practice the preconditioning matrix H is constructed to be an 
approximation to having the following properties: 

(1) H has a sparse matrix structure; 

(li) H is easily Invertible. 

One choice for H considered here is an approximate LU decomposition of a 
finite difference approximation to (2.1), l.e.. 
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H = LU 

where L and U are lower and upper triangular matrices respectively. 

Although in this section we have considered the general nonlinear equation 
(2.1) we shall restrict ourselves to the linear case in the remainder of this 
paper. Treatment of the nonlinear case will be the subject of future work. 


3. An Heuristic Discussion of Spectral Multigrid 

Of course, far better schemes are available than the simple Rlchardson"s 
method. Multigrid methods have demonstrated their ability to accelerate many 
types of relaxation schemes (see, e.g., the conference proceeedings edited by 
Hackbusch and Trottenberg (1982)). Multigrid methods have recently been 
developed for spectral discretizations of the linear version of (2.1) by Zang, 
Wong, and Hussalni (1982, 1983). Streett, Zang, and Hussaini (1983) have 

shown that these techniques are extremely effective for the nonlinear 
potential flow problem of transonic aerodynamics. Preconditioning techniques 
play a crucial role in these spectral multigrid schemes. Our purpose here is 
to investigate additional choices for the preconditioning matrix. 

The fundamentals of spectral multigrid are perhaps easiest to grasp for 
the simple model problem 



(3.1) 


on [0,2it] with periodic boundary conditions. The Fourier approximation to 
the left-hand side of (3.1) at the collocation points x^ = Zirj/N is 


N/2-1 


p=-N/2+l 


Ur, 


ipx4 

e 


(3.2) 
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where Up are the Fourier coefficients of u. The eigenfunctions of this 
approximation are 

E.(p) - , (3.3) 

with the corresponding eigenvalues 


A(P) = P^, (3.4) 

where j = and p = -N/2+1 , • • • ,N/2-l. The index p has a 

natural interpretation as the frequency of the eigenfunction. 

Consider now the iterative process described by (2.3) with H taken to be 
the identity nmitrix, i.e., without any preconditioning. The iteration matrix, 
C, of this scheme is given by 

C = I - 0) L . 

sp 


The Iterative scheme is convergent if the eigenvalues, X, of 


L satisfy 
sp 


|1 - o)X| < 1. 


The best choice of m Is that for which 


(1 - oj X ) == - (1 - to 
max 


4 ), 

min 


2 

where ^ ^min^~ largest and smallest eigenvalues 

of L„_ respectively, for then the largest values of p = 1 - toX are equal 
°P 

in magnitude and have opposite sign (see Fox (1964)). One need not worry 
about the p = 0 eigenfunction since it corresponds to the mean level of the 
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solution, which is at one's disposal for this problem. The optimal relaxation 
parameter for this single-grid procedure is 


^SG X + X . 

max min 

It produces the spectral radius 


_ max min 

’^SG X + X . 

max min 


(3.5) 


(3.6) 


2 9 

Unfortunately, ~ ^ > which implies that 0(N ) iterations are 

required to achieve convergence. 

This slow convergence is the outcome of balancing the damping of the 
lowest frequency eigenfunction with that of the highest frequency one in the 
mlnimax problem described above. The multigrid approach takes advantage of 
the fact that the low frequency modes (Ipl < N/4) can be represented just as 
well on coarser grids. It settles for balancing the middle-frequency one 
(Ipl = N/4) with the highest frequency one (Ipl = N/2), and hence damps 
effectively only those modes which cannot be resolved on coarser grids. In 
(3.5) and (3.6), is replaced by ^ A (N/4). The optimal relaxation 

parameter in this context is 


0 ). 


MG A + A , , 
max mid 


The multigrid smoothing factor 


max mid 

^MG A + A ^ 
max mid 


(3.7) 


(3.8) 
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measures the damping rate of the high-frequency modes. Alternatively, we may 
write 





1 

T 


where k:, = X /X . , is known as the multigrid condition number. In this 
example = 0.60, Independent of N. The price of this effective damping 

of the high-frequency errors is that the low-frequency errors are hardly 
damped at all. Table I compares the single-grid and multigrid damping factors 
for N = 64. However, on a grid with N/2 collocation points, the modes for 
Ipl € [N/8, N/4] are now the high-frequency ones. They get damped on this 
grid. Still coarser grids can be used until relaxations are so cheap that one 
can afford to damp all the remaining modes, or even to solve the discrete 
equations exactly. For the case illustrated in Table I the high-frequency 
error reduction in the multigrid context is roughly 250 times as fast as the 
single-grid reduction for N = 64. 


Table I. 

Dampllng Factors for 

N = 64 

1 

.9980 

.9984 

2 

.9922 

.9938 

4 

.9688 

.9750 

8 

.8751 

.9000 

12 

.7190 

.7750 

16 

.5005 

.6000 

20 

.2195 

.3750 

24 

.1239 

.1000 

28 

.5298 

.2250 

32 

.9980 

.6000 




Morcholsne (1979) and Orszag (1980) have proposed a preconditioning for 
spectral methods which amounts to using a low-order finite difference 
approximation for H. Let and L denote second-order, fourth- 

bp 

order and spectral discretizations of the operator - d /dx . The eigenvalues 
of these discretizations are given below: 


(2) _ 2[1 - cos(kAx)] 
(Ax) 


(4) _ cos(2kAx) - 16 cos(kAx) + 15 


6(Ax)‘ 


x; ^ = k . 

k 


The effective eigenvalues of the preconditioned iterations based on 
^ and ^ are then given by 




(kAx) 


2[1 - cos(kAx)] ’ 


A 


= (k^)(X 


k '' 


6(kAx)^ 

cos(2kAx) - 16 cos(kAx) + 15 
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Table II. Properties of Finite Difference Preconditioning 
for the Model Problem 


Finite Difference Order 

A . 

min 

^mld 

A 

max 

^SG 

^MG 

2 

1.00 

1.23 

2.47 

0.424 

0.336 

4 

1.00 

1.06 

1.85 

0.298 

0.273 

6 

1.00 

1.02 

1.63 

0.240 

0.231 


Similar results for even higher-order finite difference preconditionings 
are straightforward but tedious. The key properties of this class of 
preconditioning are given in Table II. It will be seen in subsequent sections 
that these results on the spread of the eigenvalues of the preconditioned 
systems agree well with those obtained computationally in two dimensions for 
Dirlchlet problems. 

2 

Unlike the original system, which has a condition number scaling as N , 
the preconditioned system has a condition number which is independent of N. 
The fourth-order finite difference operator offers around a 20% improvement in 
convergence rate over the second-order operator. This is partially offset by 
the additional cost of inverting the finite difference operator. The higher- 
order preconditionings are of doubtful utility. 

The preconditionings are clearly most effective for the longer wavelength 
eigenfunctions, as reflected by how close is to 1. In fact, for the 
fourth-order version, is already so close to 1 that the multigrid 
convergence rate is only slightly faster than the single-grid rate. The 
advantage of multigrid for these preconditioned systems only shows up in two 
dimensions. Unlike the one-dimensional case, the inversion of the two 
dimensional operator is nontrivial. It seems advisable in these situations to 
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do only an approximate Inversion of the finite difference operator, for 
example, by using an Incomplete LU decomposition of H as the actual 

preconditioner. The outcome of this choice is that while A can be kept 

max ^ 

well under control, exhibiting a very slow growth with N, plunges 

precipitously to zero. Fortunately A^^^ remains virtually unchanged . Thus 
multigrid is attractive for these two-dimensional problems . 

We describe the multigrid process by considering the interplay between two 
grids. The fine grid problem can be written in the form 


The decision to switch to the coarse grid is made after the fine grid 
£ 

approximation V has been sufficiently smoothed by the relaxation process, 

f f 

l.e., after the high-frequency content of the error V - U has been 

sufficiently reduced. The auxiliary equation on the coarse grid is 

= F^, 

where 

c f f f 

F = R[F - L V ] , 

The restriction operator R interpolates a function from the fine grid to the 
coarse grid. The coarse grid and correction are denoted by and U^, 

respectively. After an adequate approximation to the coarse grid problem 

has been obtained, the fine grid approximation is updated using 


f f c 
V -f- V + PV • 
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The prolongation operator P Interpolates a function from the coarse grid to 
the fine grid. 

The natural prolongation operator in this application represents 
trigonometric interpolation. We describe this process below. 

On the coarse grid the discrete Fourier coefficients of the corrections 
Uj at the collocation points x j are computed using 

. , -Ipx . 

"p"iT “j ' P - -N /2,. ••,1)^/2 - 1. 

^ c j=0 ■’ 


The fine’, grid approximation is then updated using 


N /2 - 1 
u . -e u , + ) 

J J p= -N^/2 


ipx , 


u e 
P 


where Xj, j=0,l,»**,N^ - 1, are the fine grid collocation points. 

The restriction operator is constructed in a similar fashion. It turns 
out that except for a factor of 2, P and R are adjoint. In the Chebyshev 
case we force R to be the adjoint of P. 


4. A Survey of Second-Order Finite Difference Preconditionings 

In this section we compare several types of preconditioning based on 
Incomplete LU decomposition (see Mei jerink and van der Vorst (1981) ) of the 
matrix which represents the standard five-point second-order finite difference 
approximation to the differential equation (2.1) . Two such preconditionings 
were discussed in Zang, Wong, and Hussalni (1983) . The first, denoted by 
has L identical to the lower triangular portion of Hpp, the finite 
difference matrix, and U chosen so that the two super diagonals of LU 
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agree with those of Hpp. The second preconditioning H^g has the diagonal 
elements of L altered from those of Hpj) so as to ensure that the row sums 
of Hpg and are identical. We introduce a third type of 
preconditioning based on the strongly implicit method of Stone (1968). 

The elements of L and U in all these instances can be easily computed 
by simple recursive formulae. The construction of the factors L and U are 
described in detail for Stone's method. 

A five-point approximation to (2.1) at the pth mesh point can be written 
as 


B 

P 


u „ + D u , 
p-N p p-1 


+ E u + F 
P P P 


>1 


+ H 

P 


'^p+N 



(4.1) 


or, in matrix form, as 

u = q . 

FD ~ ~ 

Stone's idea was to modify the matrix Hp^ by a 'small' matrix M so 
that : 

(i) the factorization of H = H + M into the product LU involves 
much less work than the standard LU decomposition of Hp^: 

(ii) the elements of L and U are easily calculated; 

(ill) IIMII « IIHppll. 

As a step to satisfying these criteria the factors L and U were chosen to 
have three non-zero diagonals, as shown in Fig. 1, corresponding to the 
diagonals B, D, and E, and E, F, and H respectively of Hpp. 

The product LU has seven non-zero diagonals, the additional two being 
immediately interior to the B and H diagonals. The matrix M is taken to 
comprise these extra two diagonals. The elements of L and U can be 
computed from their relationships with the elements of Hpj^. 
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Figure 1 


Analogous with (4.1) the pth component of Mu can be written in the form 


C u . + G u . . 
p p-N+1 p p+N-1 


(4.2) 


Stone decided to diminish the magnitude of this term by subtracting from it a 
closely equivalent expression obtained by Taylor series expansions. Using 
these expansions it is easily shown that 


u 1 ~ “U +u,, + u 

p-N+1 — p p+1 p-N 


and 


(4.3) 


u . 1 ~ “ U + u „, + U 1 . 

P+N-1 — p p+N P“1 


Stone then Introduced a parameter a, 0 < a < 1, and defined the pth component 


of Mu to be 
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C {u - af-u + u + u „!} 

pi p-N+1 p p+1 p-N^^ 

+ G {u ^ - a(-u + u + u -)}. (4.4) 

p'- p+N-1 '' p p+N p-l-'J 

Hence the pth equation of Hu = ^ is, by (4.1) and (4.4), 

(B - aC )u „ + (D - aG )u . + {e + a(C + G )}u + (F - aC )u ,, 

' p p p-N p p p-1 ‘p p P^P p p p+l 

+ (H - aG )u + C u „ T + G u = q . (4.5) 

p p p+N p p-N+1 p p+N-1 ^p ^ ' 

The relationships between the elements of H = H + M and those of L 

r u 

and U are then given by 


b = B /fl + af 
p p '■ p-N^ 

d = D /f 1 + ah , ] 

p p ^ p-l"" 




p-N 


+ d 


P P- 


h 

p-1'' 




+ d 

P 



(4.6) 


f = (f - ab f „)/e 

p ^ p p p-N^ p 

h = (h -ad h ,)/e , 
P ^ P P p-1^ P 


2 

for p = !,•••, (N-1) . Any terms with non-positive subscripts that occur in 
(4.6) are replaced by zero. For any fixed value of a, 0 < a < 1, (4.6) 
defines an incomplete LU-decompositlon of Hpp. We denote this factorization 


by Hg^(a). 
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The eigenvalues of the iteration matrices H corresponding to these 
types of preconditioning have been computed numerically by the QR algorithm 
(see Wilkinson (1965) ) . The extreme eigenvalues for Hpp, and are 
given In Zang., Wong, and Hussalni (1983) , and those for Hpp, Hgp(l.O) and 
Hgrj,(0.9) in Table III. A few of the eigenvalues at the lower end of the 
spectrum have small imaginary parts while the rest are completely real. 


Table III. Extreme Eigenvalues 

for Preconditioned Chebyshev Operator 

N 

4 

L 


.0)L 

H~J(0.9)L 


^min 

X 

max 

^min 

X 

max 

X , X 

min max 

4 

1.00 

1.76 

1.01 

1.64 

0.99 1.65 

8 

1.00 

2.13 

0.78 

2.04 

0.80 2.07 

16 

1.00 

2.31 

0.62 

2.28 

0.55 2.33 

24 

1 

1.00 

2.36 

0.58 

2.95 

0.36 2.41 



Table IV. 

Single-grid Condition Number 



-1 

-1 

-1 

-1 

N 



Hg^(1.0)L 

Hg^(0.9)L 

4 

1.85 

1.72 

1.63 

1.67 

8 

3.91 

2.71 

2.61 

2.59 

16 

11.62 

4.07 

3.67 

4.24 

24 

24.66 

5.22 

5.14 

6.69 
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Table V. Multigrid Condition Number 


N 


h;' L 

RS 

Hg^(1.0)L 

Hg^(0.9)L 

8 

1.79 

2.07 

1.70 

1.71 

16 

2.12 

2.92 

2.09 

2.08 

24 

2.26 

3.79 

2.81 

2.15 


In order to examine the effectiveness of these preconditionings from the 
multigrid point of view, we need to know the smallest high-frequency 
eigenvalue. The numerical results Indicate that this is 1.22 for and 

1.45 for H|^g, Independent of N. For Hg.p this value lies between 1.10 and 

1.21. Tables IV and V contain the single-grid and multigrid condition numbers 
respectively for %S’ Hgrj.(1.0) and Hgrj,(0.9). Here we see that 

Hg,j,(0.9) is more effective as a preconditioner for multigrid iterations. The 

multigrid condition number for Hg,p(1.0) lies between those for and 


%S* 


In Table III we see that the maximum eigenvalue of H (0.9)L grows 

Ox S p 

more slowly than that of '^ST^^ **^^^sp increasing N. For values of N 

greater than 24 the eigenvalue calculation using the QR algorithm becomes too 
expensive. Using the power method the maximum eigenvalue of Hg^(0.9)L^p 
for N = 32 was calculated and found to be 2.45. Thus the maximum eigenvalue 
of Hg^(0.9)L^^ exhibits slow growth with N. This means that the multigrid 
condition number does not increase drastically with N. However, it is 
doubtful whether the choice of a in Stone's algorithm is problem independent 


and this may detract from the robustness of this preconditioning. 
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To complete this section we look at one other type of Incomplete LU" 

decomposition which is due to Wesseling (1982a). In this decomposition the 
sparsity of the factors L and U differ from those in Fig. 1 by the 

addition of extra non-zero diagonals c and g which are located immediately 

interior to the h and h diagonals respectively. The main diagonal of U 

is specified to be unity. The elements of L and U are computed from those 
of Hpp (see (4.1)) recursively as follows: 


b 

P 


B 

P 


-b f 


p p-N 


d 

P 


D 

P 




e 

P 




*^p ®p-N+l 




(4.7) 


f = (f - h , ^ c )/e 
p ^ p p-N+1 p-* p 


g = - d h 1 /e 
P P p-1 P 


h = H /e , 
P P P 


for p = !,•••, (N-1) . Quantities that are not defined are replaced by zero. 
This preconditioning will be known as the seven-diagonal preconditioning and 
we denote it by Hgj^. The error matrix M contains two non-zero diagonals. 
It can be shown that the norm of the error matrix of the seven-diagonal 
preconditioning is smaller than those corresponding to and 
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In Table VI we give details of the seven-diagonal preconditioning for 

different values of N. Included in this table are the extreme eigenvalues of 

-1 

HgD and the single-grid and multigrid condition numbers. The smallest 

high frequency eigenvalue was found to be 1.23 Independent of N. The results 


indicate that Hgp is the most effective preconditioner considered to date. 


Table VI. Details of Seven-diagonal Preconditioning 


N 

A . 

min 

A 

max 

'^SG 

'^MG 

4 

1.00 

1.76 

1.77 

- 

8 

0.85 

2.16 

2.54 

1.76 

16 

0.46 

2.38 

5.22 

1.93 

24 

0.25 

2.47 

9.81 

2.01 


5. Higher-order Finite Difference Preconditioning 

Here we consider the possibility of choosing the preconditioning matrix 
H to be a fourth-order finite difference representation of (2.1). For a 
general non-uniform grid a compact nine-point finite difference approximation 
cannot be constructed since the set of equations for the coefficients of the 
scheme is inconsistent. Instead an approximation was constructed based on 
fourth-order finite difference formulae for each of the second derivatives 
separately. The coefficients of the function values at the grid points are 
functions of the mesh lengths which define the non-uniform mesh and are given 
in the Appendix. At internal points the approximation to the second 

derivative reduces to the following in the case when the mesh is uniform: 
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d^u _ 1 

2 2 
dx 12h 


u „ + 16u. , - 30u. + 16u. , T - u. + O(h^) . (5.1) 

1-2 1—1 1 i+i 1+2^ 


The approximation at points one mesh length from the boundary is 
constructed using the same number of points as interior equations thus 
maintaining the same order of accuracy. However , in doing this the s)raimetry 
of the sparsity pattern of the coefficient matrix is destroyed. As before, it 
is an incomplete LU decomposition of this finite difference matrix that is 
used to precondition (2.3). An algorithm due to Wessellng (1982b) was used to 
perform this decomposition. The factors L and U have the same sparsity 
pattern as the lower and upper triangular portions of Hpp respectively and 
the corresponding diagonal elements of L and U are equal. We let Hy 
denote this preconditioning. 


Table VII. Extreme Eigenvalues for Preconditioned Chebyshev Operator 


N 


^FD ^sp 


L 

sp 

A . 
min 

X 

max 

^min 

A 

max 

8 

1.00 

1.59 

0.51 

1.75 

16 

1.00 

1.79 

0.19 

2.03 

24 

1.00 

1.83 

0.09 

2.13 


Table VIII. Condition Numbers 


N Single-grid . Multigrid 


8 

3.43 

1.73 

16 

10.68 

2.01 

24 

23.67 

2.11 
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Table VII contains the largest and smallest eigenvalues of and 

-1 

^sp’ Table VIII presents the single-grid and multigrid condition numbers 


for the matrix ^sp‘ 


We note that the multigrid condition numbers 
presented here compare favourably with those obtained using second-order 
finite difference preconditioning. The smallest high frequency eigenvalue was 
found to be 1.01, Independent of the value of N. A disadvantage of this 
preconditioning Is that more work Is needed to perform the decomposition. 


6. Spectral Preconditioning 

The final type of preconditioning we investigate is that based on the 
spectral matrix, Lgp, Itself. As stated earlier the spectral matrix is full 
and hence costly to invert. Let be the matrix containing five diagonals 
of Tgp, the positions of which correspond to the non-zero diagonals of the 
finite difference matrix based on the five point formula. Let the matrix Hg, 
Illustrated in Fig. 2, be the corresponding matrix with nine non-zero 
diagonals. 

Experimentally, we found that several eigenvalues of L^^ were 
negative, thus showing that this matrix is not positive definite. The 
eigenvalues of L^p were also computed and found to be positive. 
Therefore, was discarded as a potential preconditioner and the usefulness 
of Hg examined further. 

An approximate LU-decompositlon of Hg is used to precondition (2.2), 
i.e., H is taken to be the product of a lower triangular matrix L and an 
upper triangular matrix U. We perform this decomposition according to the 
criterion described In Section 5 for constructing the factors L 
and U are chosen so that L is identical to the lower triangular portion 
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of Hg and U Is such that the four super diagonals of LU agree with those 
of Hg. The sparsity pattern of the factors L and U is shown in Fig. 3. 
We denote this preconditioning by ^gp* 

The elements of the matrices L and U can be computed recursively as 
follows: 


a = A 
P P 


b = B - a k 
p p p P-2N 


c = C 
P P 


d = D - c f . 
P P P P-1 


e = E “d f , - c g o“b k ^-a i „„ 
P P P p-1 P p-2 P p-N p p-2N 


( 6 . 1 ) 


f = (f - d g 1 1 /e 
P ^ P P P-I-^ P 


8 “ 

P P P 


k = (k - b £ J/e 

p ^ p p p-N^ p 


£ = L /e , 

P P P 


for p == 1,***,(N-1) . Any terms with non-positive subscripts occurring in 
(6.1) are replaced by zero. 


h"^ I, 
sp sp 


Table IX contains the largest and smallest eigenvalues of and 

Table X presents the single-grid and multigrid condition 


.-1 


ntimbers for the matrix ^sp’ smallest high frequency eigenvalue was 

found to be 0.90, independent of the value of N. For N = 32 we were able 

to calculate: that the maximum eigenvalue of H ^ L was 1.52. We conclude 

sp sp 

from these results that H„_ is likely to be an effective preconditioner in 

bp 
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Table IX. Extreme Eigenvalues for the Preconditioned Chebyshev Operator 


N 


L 

9 sp 


H ^ L 
sp sp 

^min 

X 

max 

^mln 

X 

max 

8 

0.63 

1.19 

0.39 

1.31 

16 

0.23 

1.33 

0.12 

1.47 

24 

0.11 

1.37 

0.05 

1.51 


Table X. Condition Numbers 


N Single-grid Multigrid 


8 

3.40 

1.44 

16 

12.53 

1.63 

24 

27.87 

1.69 


7. Numerical Results 

Here we investigate the performance of various preconditioning ideas 
developed in this paper within the SMG method. The operation count for a 
relaxation sweep of SMG is 0( N log M ) where M = (N-1)^ compared with 
0( M ) for finite difference methods. Thus we make our comparisons in terms 
of machine time instead of the work units of Brandt (1977). 

The measure used is the equivalent smoothing rate, defined by This 

was introduced in earlier work on SMG by Zang, Wong, and Hussainl (1983) and 
is defined as follows. In some preliminary calculations, the average time 
Tq required for a single fine-grid relaxation is determined. For an actual 
multigrid calculation let r^ and T 2 be the residuals after the first and 
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last fine-grid relaxations respectively, and let t be the total CPU time. 

Then we define y by 
e 

1/((t/Tq) - 1) 


There are many variants of the multigrid method. In the one used here we 
first solve the problem on the coarsest grid, that solution Is then 
interpolated to the next finer level to serve as the Initial approximation for 
a multigrid iteration involving these two levels, etc. The sizes of the grids 
on the coarsest and finest levels are 4x4 and 32 x 32 respectively. 
Internal checks based on the anticipated smoothing rates were used in 
governing decisions to switch levels. Non-stationary Richardson iteration 
employing three distinct parameters was used for relaxation. We used the 
correction scheme of Brandt (1977) with random numbers for the initial guess. 
On lower levels the right-hand sides were obtained by applying the appropriate 
restriction operator to the finest level right-hand side. 

The test problems are specified by 



a(x,y) = b(x,y) = 1 + e exp[ cos(gTr(x+y))) , 
u(x,y) = sln(aTrx + it/4) sin(aiTy + tt/4). 

The parameters of the test problem are given in Table XI. Problem 1 has 
constant coefficients and is also well-resolved by the Chebyshev collocation 
method. On the other hand, for Problem 3 the coefficients of the equation 
oscillate so rapidly that the finest grid cannot resolve them. The converged 
solution of the finest collocation equations has an error of order 1. This is 
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a test of whether the Chebyshev SMG method is robust enough to converge on 
such a problem. Table XII contains the values of the equivalent smoothing 
rates for various preconditioners. 


Table XI. Parameters of the Test Problems 


Problem 

e 

a 


1 

0.00 

1 

1 

2 

0.20 

2 

2 

3 

1.00 

5 

10 


Table 

XII. 

Equivalent Smoothing Kates 


Problem 

%u 

Hg^(0.9) 

^SD 

Hsp 

1 

0.26 

0.16 

0.12 

0.24 

2 

0.60 

0.65 

0.47 

0.54 

3 

0.76 

0.84 

0.76 

0.82 


The results show the seven-diagonal finite difference preconditioning and 
the spectral preconditioning are comparable in terms of efficiency and give an 
improvement over the other preconditioners except on Problem 3 where the 
performance of all the preconditioners is about the same. However, since in 
our computations to date the spectral matrix is never stored, the elements 
of Hgp are more expensive to compute than Notice the superb 
performance of Stone's preconditioning with a = 0.9 on the constant 
coefficient problem. This performance is not maintained on the more difficult 
problems which demonstrates that the choice of the parameter a is problem 
dependent. The seven-diagonal preconditioning also performs extremely well on 
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the constant coefficient problem. We did not experiment with the fourth-order 
finite difference preconditioning since to compute the approximation in the 
non-constant coefficient case was thought to be too laborious. 

8 . Conclusions 

In this paper we have developed efficient iterative techniques for solving 
the algebraic equations which arise from the discretization of a self-adjoint 
elliptic equation using the spectral collocation method. An advantage of 
using spectral methods is that they possess superior approximation properties 
compared with finite difference and finite element methods. In practice, this 
means one can obtain the same accuracy with fewer mesh points than are needed 
for finite difference or finite element methods. Acceleration of the basic 
iterative scheme has been enhanced by employing the multigrid method. The 
need for preconditioning has been demonstrated and efficient preconditioners 
for multigrid Iterations presented. 
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Appendix 

Higher-order Finite Differences on Non-uniform Grids 

Suppose that we require an approximation to the second derivative of u 
at point C in Fig. 4. We assume that the approximation has the form 

d^u 

dx 

where the coefficients a, b, c, d and e are to be determined. 

p q r s 

I I I I I 

ABC D E 


Figure 4 

The right-hand side of (A.l) is expanded in a Taylor's series about the 
point C. The coefficient of the second derivative is set to unity while 
those of u and the first, third, and fourth derivatives are set to zero. 
This results in a system of five linear equations for the coefficients 

O 

appearing in (A.l). The accuracy of the approximation is 0(h) where h is 
a typical mesh length. 

If we define G = p + q and H = r + s, then the coefficients are given 


by 
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-2q(H + r) + 2Hr 
[Gp(G + H)(r + G) ] ’ 


2G(H + r) - 2Hr 
[pq(q + r)(H + q)] ’ 


2H(G + q) - 2Gq 
[rs(q + r)(G + r) ] ’ 


-2r(G + q) + 2Gq 
[Hs(H + q)(H + G)] ’ 
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